
# load ipe package (REQUIRES 64-BIT WINDOWS DUE TO C++ CODE COMIPLED AS A DLL)
install.packages("https://sites.google.com/a/stonybrook.edu/mperess/ipe_PA.tar.gz")

# some helper functions
source("r0.r")

# load the needed libraries
library0("ipe")
library0("wnominate")
library0("pscl") # for ideal
library0("emIRT")

# replicate the monte carlos
set.seed(1234)
Ns <- c(200,500,1000,2000)
Ts <- c(200,500,1000,2000)
cols <- c('blue','red','magenta','cyan')

Alpha0Mat <- list()
AlphaMat <- list()
AlphaSeMat <- list()
for(i in 1:length(Ns))
{
	N <- Ns[i]
	T <- Ts[i]
	Alpha0 <- -2 + 4 * runif(N)
	Delta00 <- rnorm(T)
	Delta10 <- 0.1 + runif(T)
	ipe0 <- list(Alpha=matrix(Alpha0,ncol=1),Delta=cbind(Delta00,Delta10))
	Party <- rep("R",N)
	Party[Alpha0 < 0] <- "D"
	ipe0 <- ipe_normalize_by_group(ipe0,Party,list(D=-1,R=1))
	ipe_means_by_group(ipe0,Party)
	Alpha0 <- ipe0$Alpha[,1]
	Delta00 <- ipe0$Delta[,1]
	Delta10 <- ipe0$Delta[,2]

	R <- 100
	AlphaMat[[i]] <- matrix(rep(NA,N*R),ncol=N)
	AlphaSeMat[[i]] <- matrix(rep(NA,N*R),ncol=N)
	Alpha0Mat[[i]] <- Alpha0

	for(r in 1:R)
	{
		print0("r",r)
		YStar <- rep(1,N) %*% t(Delta00) + Alpha0 %*% t(Delta10) + matrix(rnorm(N*T),nrow=N)
		Miss <- matrix(runif(N*T),nrow=N)<.7 # ensures ~70% missing
		Y1 <- ((YStar >= 0) + 1) * !Miss 
		Y1s <- ipe_dense_to_sparse(Y1)

		ipe1s <- ipe_start_sparse(Y1s,1)
		ipe1u <- ipe_pmle_sparse(Y1s,1,ipe=ipe1s,InferenceOpt=1)
		ipe1  <- ipe_normalize_by_group(ipe1u,Party,list(D=-1,R=1))

		AlphaMat[[i]][r,] <- ipe1$Alpha
		AlphaSeMat[[i]][r,] <- ipe1$AlphaSe

		gc()
	}
}

############
# FIGURE 1 #
############

pdf(file="fig1.pdf",width=6,height=6)
plot(0:1,type='n',xlim=c(-2,2),ylim=c(-1,1),xlab="Alpha",ylab="Bias")
for(i in 1:4)
{
	Bias <- colMeans(AlphaMat[[i]])-Alpha0Mat[[i]]	
	points(lowess(Bias ~ Alpha0Mat[[i]]),type='l',col=cols[i])
}
legend(.75,1,c("N=T=200","N=T=500","N=T=1000","N=T=2000"),lty=c(1,1,1,1),col=cols)
dev.off()

############
# FIGURE 2 #
############

pdf(file="fig2.pdf",width=6,height=6)
plot(0:1,type='n',xlim=c(-2,2),ylim=c(0,.5),xlab="Alpha",ylab="RMSE")
for(i in 1:4)
{
	RMSE <- colMeans((AlphaMat[[i]] - matrix(rep(1,R),ncol=1) %*% Alpha0Mat[[i]])^2)^.5
	points(lowess(RMSE ~ Alpha0Mat[[i]]),type='l',col=cols[i])
}
legend(.75,.5,c("N=T=200","N=T=500","N=T=1000","N=T=2000"),lty=c(1,1,1,1),col=cols)
dev.off()

############
# FIGURE 3 #
############

pdf(file="fig3.pdf",width=6,height=6)
plot(0:1,type='n',xlim=c(-2,2),ylim=c(.85,1.05),xlab="Alpha",ylab="95% Coverage")
for(i in 1:4)
{
	Low <- AlphaMat[[i]] - 1.96 * AlphaSeMat[[i]]
	High <- AlphaMat[[i]] + 1.96 * AlphaSeMat[[i]]
	Cov <- colMeans(Low <= matrix(rep(1,R),ncol=1) %*% Alpha0Mat[[i]] & matrix(rep(1,R),ncol=1) %*% Alpha0Mat[[i]] <= High)
	points(lowess(Cov ~ Alpha0Mat[[i]]),type='l',col=cols[i])
}
abline(h=.95,lty=2)
legend(.75,1.05,c("N=T=200","N=T=500","N=T=1000","N=T=2000"),lty=c(1,1,1,1),col=cols)
dev.off()

# 112th Senate
rc <- readKH("sen112kh.ord")
Y <- ipe_rc_to_dense(rc)
party <- rc$legis.data$party

# w-nominate
ptm <- proc.time()
wnom <- wnominate(rc,dims=1,polarity=c(2)) # polarity based on Sessions
tm11 <- proc.time() - ptm
tm11

# w-nominate, bootstrap
ptm <- proc.time()
wnom <- wnominate(rc,trials=20,dims=1,polarity=c(2)) # polarity based on Sessions
tm12 <- proc.time() - ptm
tm12
ipe_wnom <- wnom_to_ipe(wnom)
ipe_wnom <- ipe_normalize_by_group(ipe_wnom,party,list(D=-1,R=1),num=30)
alpha_wnom <- ipe_wnom$Alpha

# ideal
ptm <- proc.time()
ideal <- ideal(rc,d=1,normalize=TRUE,maxiter=5000,burnin=2500,thin=1)
tm13 <- proc.time() - ptm
tm13
ipe_ideal <- ideal_to_ipe(ideal,rc)
ipe_ideal <- ipe_normalize_by_group(ipe_ideal,party,list(D=-1,R=1),num=30)
alpha_ideal <- ipe_ideal$Alpha

# using PMLE, no inference
ptm <- proc.time()
ipe1s <- ipe_start_dense(Y,D=1)
ipe1u <- ipe_pmle_dense(Y,D=1,ipe=ipe1s,SkipInference=1)
tm14 <- proc.time() - ptm
tm14

# using PMLE, inference
ptm <- proc.time()
ipe1s <- ipe_start_dense(Y,D=1)
ipe1u <- ipe_pmle_dense(Y,D=1,ipe=ipe1s)
ipe1 <- ipe_normalize_by_group(ipe1u,party,list(D=-1,R=1),num=30)
ipe1$Alpha[ipe1$IndFit[,1]<30] <- NA
tm15 <- proc.time() - ptm
tm15
alpha_pmle <- ipe1$Alpha

# using emIRT
ptm <- proc.time()
rc2 <- convertRC(rc)
p <- makePriors(rc2$n, rc2$m, 1)
s <- getStarts(rc2$n, rc2$m, 1)
emirt <- binIRT(.rc = rc2,.starts = s,.priors = p,.control = {list(threads = 1,verbose = FALSE,thresh = 1e-6)})
tm16 <- proc.time() - ptm
tm16

# using emIRT, bootstrap
ptm <- proc.time()
p <- makePriors(rc2$n, rc2$m, 1)
s <- getStarts(rc2$n, rc2$m, 1)
emirt <- binIRT(.rc = rc2,.starts = s,.priors = p,.control = {list(threads = 1,verbose = FALSE,thresh = 1e-6)})
emirt <- boot_emIRT(emirt, .data = rc2, .starts = s, .priors = p,.control = list(threads = 1, verbose = FALSE, thresh = 1e-06), Ntrials=20, verbose=2)
tm17 <- proc.time() - ptm
tm17
ipe_emirt <- emirt_to_ipe(emirt,rc2)
ipe_emirt <- ipe_normalize_by_group(ipe_emirt,party,list(D=-1,R=1),num=30)
alpha_emirt <- ipe_emirt$Alpha

############
# Figure 4 #
############

pdf(file="fig4.pdf",width=6,height=6)
pairs(~alpha_wnom+alpha_ideal+alpha_emirt+alpha_pmle,labels=c("wnominate","ideal","emIRT","ipe"),main="",xlim=c(-1,2),ylim=c(-1,2))
dev.off()

# 112th House
rc <- readKH("hou112kh.ord")
Y <- ipe_rc_to_dense(rc)
party <- rc$legis.data$party

# w-nominate
ptm <- proc.time()
wnom <- wnominate(rc,dims=1,polarity=c(2)) # polarity based on Sessions
tm21 <- proc.time() - ptm
tm21

# w-nominate, bootstrap
ptm <- proc.time()
wnom <- wnominate(rc,trials=20,dims=1,polarity=c(2)) # polarity based on Sessions
tm22 <- proc.time() - ptm
tm22

# ideal
ptm <- proc.time()
ideal <- ideal(rc,d=1,normalize=TRUE,maxiter=5000,burnin=2500,thin=1)
tm23 <- proc.time() - ptm
tm23

# using PMLE, no inference
ptm <- proc.time()
ipe1s <- ipe_start_dense(Y,D=1)
ipe1u <- ipe_pmle_dense(Y,D=1,ipe=ipe1s,SkipInference=1)
tm24 <- proc.time() - ptm
tm24

# using PMLE, inference
ptm <- proc.time()
ipe1s <- ipe_start_dense(Y,D=1)
ipe1u <- ipe_pmle_dense(Y,D=1,ipe=ipe1s)
ipe1 <- ipe_normalize_by_group(ipe1u,party,list(D=-1,R=1),num=30)
ipe1$Alpha[ipe1$IndFit[,1]<30] <- NA
tm25 <- proc.time() - ptm
tm25

# using emIRT
ptm <- proc.time()
rc2 <- convertRC(rc)
p <- makePriors(rc2$n, rc2$m, 1)
s <- getStarts(rc2$n, rc2$m, 1)
emirt <- binIRT(.rc = rc2,.starts = s,.priors = p,.control = {list(threads = 1,verbose = FALSE,thresh = 1e-6)})
tm26 <- proc.time() - ptm
tm26

# using emIRT, bootstrap
ptm <- proc.time()
p <- makePriors(rc2$n, rc2$m, 1)
s <- getStarts(rc2$n, rc2$m, 1)
emirt <- binIRT(.rc = rc2,.starts = s,.priors = p,.control = {list(threads = 1,verbose = FALSE,thresh = 1e-6)})
emirt <- boot_emIRT(emirt, .data = rc2, .starts = s, .priors = p,.control = list(threads = 1, verbose = FALSE, thresh = 1e-06), Ntrials=20, verbose=2)
tm27 <- proc.time() - ptm
tm27

# 90th Senate
rc <- readKH("sen90kh.ord")
Y <- ipe_rc_to_dense(rc)
party <- as.character(rc$legis.data$party)
state <- as.character(rc$legis.data$state)
south <- state == "AL" |  state == "AR" |  state == "FL" |  state == "GA" |  state == "LA" |  state == "MS" |  state == "NC" |  state == "SC" |  state == "TN" |  state == "TX" |  state == "VA" |  state == "WV"
party2 <- party
for(i in 1:length(party2)) if(party2[i] == "D") if(south[i] == 1) party2[i] <- "SD"

# w-nominate
ptm <- proc.time()
wnom <- wnominate(rc,dims=2,polarity=c(2,2)) # polarity based on Sessions
tm31 <- proc.time() - ptm
tm31

# w-nominate, bootstrap
ptm <- proc.time()
wnom <- wnominate(rc,trials=20,dims=2,polarity=c(2,2)) # polarity based on Sessions
tm32 <- proc.time() - ptm
tm32
ipe_wnom <- wnom_to_ipe(wnom)
ipe_wnom <- ipe_normalize_by_group(ipe_wnom,party2,list(D=c(-1,-.25),SD=c(-.6,.25),R=c(1,0)),num=30)
alpha_wnom <- ipe_wnom$Alpha

# ideal
ptm <- proc.time()
ideal <- ideal(rc,d=2,normalize=TRUE,maxiter=5000,burnin=2500,thin=1)
tm33 <- proc.time() - ptm
tm33
ipe_ideal <- ideal_to_ipe(ideal,rc)
ipe_ideal <- ipe_normalize_by_group(ipe_ideal,party2,list(D=c(-1,-.25),SD=c(-.6,.25),R=c(1,0)),num=30)
alpha_ideal <- ipe_ideal$Alpha

# using PMLE, no inference
ptm <- proc.time()
ipe1s <- ipe_start_dense(Y,D=2)
ipe1u <- ipe_pmle_dense(Y,D=2,ipe=ipe1s,SkipInference=1)
tm34 <- proc.time() - ptm
tm34

# using PMLE, inference
ptm <- proc.time()
ipe1s <- ipe_start_dense(Y,D=2)
ipe1u <- ipe_pmle_dense(Y,D=2,ipe=ipe1s)
ipe1 <- ipe_normalize_by_group(ipe1u,party2,list(D=c(-1,-.25),SD=c(-.6,.25),R=c(1,0)),num=30)
ipe1$Alpha[ipe1$IndFit[,1]<30,] <- NA
tm35 <- proc.time() - ptm
tm35
alpha_pmle <- ipe1$Alpha

# emIRT does not support two dimensions
tm36 <- rep(NA,3)
tm37 <- rep(NA,3)

############
# FIGURE 5 #
############

pdf(file="fig5a.pdf",width=6,height=6)
pairs(~alpha_wnom[,1]+alpha_ideal[,1]+alpha_pmle[,1],main="First Dimension",labels=c("wnominate","ideal","ipe"),xlim=c(-2,2),ylim=c(-2,2))
dev.off()

pdf(file="fig5b.pdf",width=6,height=6)
pairs(~alpha_wnom[,2]+alpha_ideal[,2]+alpha_pmle[,2],main="Second Dimension",labels=c("wnominate","ideal","ipe"),xlim=c(-.7,.7),ylim=c(-.7,.7))
dev.off()

###########
# TABLE 1 #
###########

tab1 <- matrix(c(tm11[3],tm12[3],tm13[3],tm13[3],tm16[3],tm17[3],tm14[3],tm15[3],tm21[3],tm22[3],tm23[3],tm23[3],tm26[3],tm27[3],tm24[3],tm25[3],tm31[3],tm32[3],tm33[3],tm33[3],tm36[3],tm37[3],tm34[3],tm35[3]),ncol=3)
colnames(tab1) <- c("112th Senate (1D)","112th House (1D)","90th Senate (2D)")
rownames(tab1) <- c("wnom","wnom w/ s.e.","ideal","ideal w/ s.e.","emIRT","emIRT w/ s.e.","ipe","ipe w/ s.e.")
tab1

# read in all the house and senate data and store in list
chambers <- list()
curr <- 1
for(j in 1:2) for(cong in 1:113)
{
	print(paste(j,cong,sep="_"))
	flush.console()
	if(j == 1)
	{
		if(cong < 10) {
			file1 <- readKH(paste("hou0",cong,"kh.ord",sep=""))
		} else {
			file1 <- readKH(paste("hou",cong,"kh.ord",sep=""))
		}
	} else {
		if(cong < 10) {
			file1 <- readKH(paste("sen0",cong,"kh.ord",sep=""))
		} else {
			file1 <- readKH(paste("sen",cong,"kh.ord",sep=""))
		}
	}
	votes <- file1$votes
	Ncurr <- nrow(votes)
	Tcurr <- ncol(votes)
	Ycurr <- 1 * (votes == 6) + 2 * (votes == 1)
	Nidcurr <- file1$legis.data$icpsrLegis
	Tidcurr <- paste(cong,"_",j,"_",1:Tcurr,sep="")
	rownames(Ycurr) <- Nidcurr
	colnames(Ycurr) <- Tidcurr
	if(j == 1) {
		inddat <- as.data.frame(cbind(rownames(file1$legis.data),as.character(file1$legis.data$state),as.character(file1$legis.data$party),rep("H",Ncurr)),stringsAsFactors=FALSE)
	} else {
		inddat <- as.data.frame(cbind(rownames(file1$legis.data),as.character(file1$legis.data$state),as.character(file1$legis.data$party),rep("S",Ncurr)),stringsAsFactors=FALSE)
	}
	colnames(inddat) <- c("name","state","party","chamber")
	inddat$party[is.na(inddat$party)] <- "NA"
	if(j == 1) {
		chambers[[curr]] <- list(Y=Ycurr,inddat=inddat,name=paste("H",cong,sep=""))
	} else {
		chambers[[curr]] <- list(Y=Ycurr,inddat=inddat,name=paste("S",cong,sep=""))
	}
	curr <- curr + 1
}

# just the senate data
chambers2 <- list()
for(i in 114:226) chambers2[[i-113]] <- chambers[[i]]

# dynamic senate, D=1

# time, no inference
ptm <- proc.time()
set.seed(1234)
ipe1s <- ipe_bridge_start(chambers2,D=1)
ipe1u <- ipe_pmle_sparse(ipe1s$Y,D=1,ipe=ipe1s$ipe,SkipInference=1)
tm41 <- proc.time() - ptm
tm41

# time, inference
ptm <- proc.time()
set.seed(1234)
ipe1s <- ipe_bridge_start(chambers2,D=1)
ipe1u <- ipe_pmle_sparse(ipe1s$Y,D=1,ipe=ipe1s$ipe)
ipe1 <- ipe_normalize_by_group(ipe1u,ipe1s$inddat$party,list("D"=-1,"R"=1),num=30)
ipe1$Alpha[ipe1$IndFit[,1]<30] <- NA
tm42 <- proc.time() - ptm
tm42

# DW-nominate
dwdat <- read_excel("dwcs.xlsx")
dwid <- dwdat$ICPSR
dw <- dwdat$Alpha1
alpha_dw <- dw[fmatch(ipe1s$Y$Nids,dwid)]
dw1 <- list()
dw1$Alpha <- matrix(alpha_dw,ncol=1)
dw1 <- ipe_normalize_by_group(dw1,ipe1s$inddat$party,list("D"=-1,"R"=1),num=0)

# emIRT
set.seed(1234)
Yd <- ipe_sparse_to_dense(ipe1s$Y)
rc <- list()
rc$n <- nrow(Yd)
rc$m <- ncol(Yd)
rc$votes <- (Yd != 0) * 2 * (Yd - 1.5)
ptm <- proc.time()
p <- makePriors(rc$n, rc$m, 1)
s <- getStarts(rc$n, rc$m, 1)
emirt1 <- binIRT(.rc = rc,.starts = s,.priors = p)
ipe3 <- ipe_normalize_by_group(emirt1$means$x,ipe1s$inddat$party,list("D"=-1,"R"=1))
tm43 <- proc.time() - ptm
tm43

# scale coordinate of emIRT for pretier plots
ipe3 <- ipe3 / 4

############
# Figure 6 #
############

pdf(file="fig6.pdf",width=9,height=5)
par(mfrow=c(1,2))
plot(dw1$Alpha,ipe3,main="emIRT vs. DWCS",xlab="DWCS",ylab="emIRT",ylim=c(-4,4))
plot(dw1$Alpha,ipe1$Alpha,main="ipe vs. DWCS",xlab="DWCS",ylab="ipe",ylim=c(-4,4))
dev.off()

# get party means in each congress
partys <- c("D","Democratic-Republican","Federalist","Jackson","R","Whig")
ipmat1 <- matrix(rep(NA,6*113),ncol=113)
ipmat2 <- matrix(rep(NA,6*113),ncol=113)
ipmat3 <- matrix(rep(NA,6*113),ncol=113)
for(i in 1:113)
{
	idx <- fmatch(rownames(chambers2[[i]]$Y),ipe1s$Y$Nids)
	for(j in 1:length(partys))
	{
		ipmat1[j,i] <- mean0(ipe1$Alpha[idx],subset=chambers2[[i]]$inddat[,3]==partys[j])
		n1 <- n0(ipe1$Alpha[idx],subset=chambers2[[i]]$inddat[,3]==partys[j])
		if(n1 < 10) ipmat1[j,i] <- NA

		ipmat2[j,i] <- mean0(dw1$Alpha[idx],subset=chambers2[[i]]$inddat[,3]==partys[j])
		n1 <- n0(dw1$Alpha[idx],subset=chambers2[[i]]$inddat[,3]==partys[j])
		if(n1 < 10) ipmat2[j,i] <- NA

		ipmat3[j,i] <- mean0(ipe3[idx],subset=chambers2[[i]]$inddat[,3]==partys[j])
		n1 <- n0(ipe3[idx],subset=chambers2[[i]]$inddat[,3]==partys[j])
		if(n1 < 10) ipmat3[j,i] <- NA
	}
}

# correct Poole and Rosenthal's coding of Democratic-Republicans
for(i in 1:30) {
	ipmat1[2,i] <- ipmat1[5,i]
	ipmat1[5,i] <- NA
	ipmat2[2,i] <- ipmat2[5,i]
	ipmat2[5,i] <- NA
	ipmat3[2,i] <- ipmat3[5,i]
	ipmat3[5,i] <- NA
}

############
# Figure 7 #
############

pdf(file="fig7.pdf",width=9,height=4)
par(mfrow=c(1,3))
cols <- c('blue','orange','green','magenta','red','cyan')
partys[1] <- "Democratic"
partys[5] <- "Republican"

plot(0:1,type='n',xlim=c(0,114),ylim=c(-2.5,4.5),xlab="Congress",ylab="",main="DW-Nominate Common Space")
for(j in 1:length(partys)) points(1:113,ipmat2[j,],col=cols[j],type='l')

legend(20,4.5,partys,lty=rep(1,6),col=cols)

plot(0:1,type='n',xlim=c(0,114),ylim=c(-2.5,4.5),xlab="Congress",ylab="",main="emIRT")
for(j in 1:length(partys)) points(1:113,ipmat3[j,],col=cols[j],type='l')

plot(0:1,type='n',xlim=c(0,114),ylim=c(-2.5,4.5),xlab="Congress",ylab="",main="ipe")
for(j in 1:length(partys)) points(1:113,ipmat1[j,],col=cols[j],type='l')
dev.off()

# dynamic sen/house, 1D, no inference #
ptm <- proc.time()
set.seed(1234)
ipe1s <- ipe_bridge_start(chambers,D=1)
ipe1u <- ipe_pmle_sparse(ipe1s$Y,D=1,ipe=ipe1s$ipe,SkipInference=1)
tm51 <- proc.time() - ptm
tm51

# dynamic sen/house, 1D, inference #
ptm <- proc.time()
set.seed(1234)
ipe1s <- ipe_bridge_start(chambers,D=1)
ipe1u <- ipe_pmle_sparse(ipe1s$Y,D=1,ipe=ipe1s$ipe)
tm52 <- proc.time() - ptm
tm52

# dynamic sen/house, 2D, no inference
ptm <- proc.time()
set.seed(1234)
ipe1s <- ipe_bridge_start(chambers,D=2)
ipe1u <- ipe_pmle_sparse(ipe1s$Y,D=2,ipe=ipe1s$ipe,SkipInference=1,MaxIter=10000)
tm61 <- proc.time() - ptm
tm61

# dynamic sen/house, 2D, inference
ptm <- proc.time()
set.seed(1234)
ipe1s <- ipe_bridge_start(chambers,D=2)
ipe1u <- ipe_pmle_sparse(ipe1s$Y,D=2,ipe=ipe1s$ipe,MaxIter=10000)
tm62 <- proc.time() - ptm
tm62

# big example
merge1 <- list()
merge1$Y <- ipe_read_sparse("merge_Y.dat")
party <- str_split_vec(merge1$Y$Nids,"_")[,2]
chamber <- str_split_vec(merge1$Y$Nids,"_")[,3]
merge1$inddat <- as.data.frame(cbind(party,chamber),stringsAsFactors=FALSE)
merge1$source <- readLines("merge_source.dat")

# split the chambers (necesary for starting values calculation)
chambers3 <- ipe_split_chamber(merge1$Y,merge1$source,inddat=merge1$inddat)

# put the chambers back together
merge2 <- ipe_merge_chambers(chambers3)
AlphaFac <- rep(1,merge2$Y$N)
AlphaFac[merge2$inddat$chamber=="V"] <- .01

# pmle, no inference
ptm <- proc.time()
set.seed(1234)
ipe1s <- ipe_bridge_start(chambers3,1,AlphaFac=AlphaFac,MinNumAlpha=10)
ipe1u <- ipe_pmle_sparse(ipe1s$Y,1,ipe=ipe1s$ipe,AlphaFac=AlphaFac,IndWeights=AlphaFac,LineSearchOpt=2,SkipInference=1)
tm71 <- proc.time() - ptm
tm71

# pmle, inference
ptm <- proc.time()
set.seed(1234)
ipe1s <- ipe_bridge_start(chambers3,1,AlphaFac=AlphaFac,MinNumAlpha=10)
ipe1u <- ipe_pmle_sparse(ipe1s$Y,1,ipe=ipe1s$ipe,AlphaFac=AlphaFac,IndWeights=AlphaFac,LineSearchOpt=2)
ipe1  <- ipe_normalize_by_group(ipe1u,ipe1s$inddat$party,list("D"=-1,"R"=1),num=10)
ipe1$Alpha[ipe1$IndFit[,1]<10] <- NA
tm72 <- proc.time() - ptm
tm72

############
# Figure 8 #
############

pdf(file="fig8.pdf",width=6,height=6)
ipe_density_by_group(ipe1,ipe1s$inddat$party,groups=c("G","D","d","i","r","R","L"),cols=c('green','blue','blue','black','red','red','yellow'),lty=c(1,1,2,2,2,1,1))
dev.off()

###########
# Table 2 #
###########

tm44 <- rep(NA,3)
tm53 <- rep(NA,3)
tm54 <- rep(NA,3)
tm63 <- rep(NA,3)
tm64 <- rep(NA,3)
tm73 <- rep(NA,3)
tm74 <- rep(NA,3)

tab2 <- matrix(c(tm43[3],tm44[3],tm41[3],tm42[3],tm53[3],tm54[3],tm51[3],tm52[3],tm63[3],tm64[3],tm61[3],tm62[3],tm73[3],tm74[3],tm71[3],tm72[3]),ncol=4)
colnames(tab2) <- c("Dynamic Senate (1D)","Dynamic House and Senate (1D)","Dynamic House and Senate (2D)","NPAT/CCES/Roll Call")
rownames(tab2) <- c("emIRT","emIRT w/ s.e.","ipe","ipe w/ s.e.")

# the tables:
write.table(tab1,file="Table1.csv",sep=",",row.names=FALSE,col.names=FALSE)
write.table(tab2,file="Table2.csv",sep=",",row.names=FALSE,col.names=FALSE)

